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Abstract — In this work, we present a novel construction 
for solving the linear multiuser detection problem using the 
Gaussian Belief Propagation algorithm. Our algorithm yields 
an efficient, iterative and distributed implementation of the 
MMSE detector. Compared to our previous formulation, the 
new algorithm offers a reduction in memory requirements, the 
number of computational steps, and the number of messages 
passed. We prove that a detection method recently proposed by 
Montanari et ai. is an instance of ours, and we provide new 
convergence results applicable to both. 

I. Introduction 

Belief propagation (BP), also known as the sum-product 
algorithm, is a powerful and efficient tool in solving, exactly 
or approximately, inference problems in probabilistic graphical 
models. The underlying essence of estimation theory is to 
detect a hidden input to a channel from its observed output. 
The channel can be represented as a certain graphical model, 
while the detection of the channel input is equivalent to 
performing inference in the corresponding graph. 

The use of BP [1] for detection purposes has been proven to 
be very beneficial in several applications in communications. 
For randomly-spread code-division multiple-access (CDMA) 
in the large-system limit, Kabashima has introduced a tractable 
BP-based multiuser detection (MUD) scheme, which exhibits 
near-optimal error performance for binary-input additive white 
Gaussian noise (BI-AWGN) channels [2]. This message- 
passing scheme has recently been extended to the case where 
the ambient noise level is unknown [3], [4]. As for sub-optimal 
detection, the nonlinear soft parallel interference cancellation 
(PIC) detector was reformulated by Tanaka and Okada as an 
approximate BP solution [5] to the MUD problem. 

In contrast to the dense, fully-connected nature of the 
graphical model of the non-orthogonal CDMA channel, a one- 
dimensional (1-D) intersymbol interference (ISI) channel can 
be interpreted as a cycle-free tree graph [6]. Thus, detection 
in 1-D ISI channels (termed equalization) can be performed 
in an optimal maximum a-posteriori (MAP) manner via BP, 
also known in this context as the forward/backward, or BCJR, 
algorithm [7]. Also, Kurkoski etal. [8], [9] have proposed an 
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iterative BP-like detection algorithm for 1-D ISI channels that 
uses a parallel message-passing schedule and achieves near- 
optimal performance. 

For the intermediate regime of non-dense graphs but 
with many relatively short loops, extensions of BP to two- 
dimensional ISI channels have been considered by Marrow 
and Wolf [10], and recently Shental etal. [11]-[13] have 
demonstrated the near-optimality of a generalized version 
of BP for such channels. Recently, BP has been proved 
to asymptotically achieve optimal MAP detection for sparse 
linear systems with Gaussian noise [14], [15], for example, in 
CDMA with sparse spreading codes. 

An important class of practical sub-optimal detectors is 
based on linear detection. This class includes, for instance, 
the conventional single-user matched filter (MF), decorre- 
lator (a.k.a. zero-forcing equalizer), linear minimum mean- 
square error (MMSE) detector and many other detectors with 
widespread applicability [16], [17]. In general, Unear detection 
can be viewed as the solution to a (deterministic) set of 
linear equations describing the original (probabilistic) estima- 
tion problem. Note that the mathematical operation behind 
linear detection extends to other tasks in communication, e.g. , 
channel precoding at the transmitter [18]. 

Recently, linear detection has been explicitly linked to 
BP [19], using a Gaussian belief propagation (GaBP) algo- 
rithm. This allows for a distributed implementation of the 
linear detector [20], circumventing the need of, potentially 
cumbersome, direct matrix inversion (via, e.g. , Gaussian elim- 
ination). The derived iterative framework was compared quan- 
titatively with 'classical' iterative methods for solving systems 
of linear equations, such as those investigated in the context 
of linear implementation of CDMA demodulation [21]-[23]. 
GaBP is shown to yield faster convergence than these standard 
methods. Another important work is the BP-based MUD, 
recently derived and analyzed by Montanari et al. [24] for 
Gaussian input symbols. 

There are several drawbacks to the linear detection tech- 
nique of [19]. First, the input matrix R„xn = ^nxk^k>:n (the 
chip correlation matrix) needs to be computed prior to running 
the algorithm. This computation requires n^fc operations. In 
case where the matrix S is sparse [15], the matrix R might 



not no longer be sparse. Second, GaBP uses 2n^ memory to 
store the messages. For a large n this could be prohibitive. 

In this paper, we propose a new construction that addresses 
those two drawbacks. In our improved construction, given a 
non-rectangular CDMA matrix S„x/£, we compute the MMSE 
detector x = (S^S + ^y^S'^y where ^I^ is the AWGN 
diagonal covariance matrix. We utilize the GaBP algorithm 
which is an efficient iterative distributed algorithm. The new 
construction uses only 2nk memory for storing the messages. 
When k <^ n this represents significant saving relative to the 
in our previously proposed algorithm. Furthermore, we do 
not explicitly compute S^S, saving an extra n^k overhead. 

We show that Montanari's algorithm [24] is an instance of 
our method. By showing this, we are able to prove new con- 
vergence results for Montanari's algorithm. Montanari proves 
that his method converges on normalized random-spreading 
CDMA sequences, assuming Gaussian signaling. Using binary 
signaling, he conjectures convergence to the large system limit. 
Here, we extend Montanari's result, to show that his algorithm 
converges also for non-random CDMA sequences when binary 
signaling is used, under weaker conditions. Another advantage 
of our work is that we allow different noise levels per bit 
transmitted. 

The paper is organized as follows. Section HIl formulates the 
problem of linear detection and presents the distributed GaBP- 
based linear detection scheme. Section |III] describes a novel 
construction for efficiently computing the MMSE detector. 
The relation to a factor graph construction is explored in 
Section |IV] New convergence results for Montanari's work 
are presented in Section |V] We conclude in Section |VI] In the 
Appendix we further explore the relation to Montanari's work. 

We shall use the following notations. The operator j }^ 
stands for a vector or matrix transpose, the matrix Ijv is a 
NxN identity matrix, while the symbols {-ji and {-Jij denote 
entries of a vector and matrix, respectively. N(j) is the set of 
graph node connected to node i. 

II. Linear Detection via Belief Propagation 

Consider a discrete-time channel with a real input vec- 
tor X. = {xi, . . . ,xk}^ governed by an arbitrary prior 
distribution, P^, and a corresponding real output vector 
y = {yi, . . . , vk}^ = /{x^} e R^. Here, the function /{•} 
denotes the channel transformation. By definition, linear de- 
tection compels the decision rule to be 

X- A{x*}- A{A-ib}, (1) 

where b = y is the K x 1 observation vector and the 
K X K matrix A is a positive-definite symmetric matrix 
approximating the channel transformation. The vector x* is 
the solution (over R) to Ax = b. Estimation is completed 
by adjusting the (inverse) matrix-vector product to the input 
alphabet, dictated by P^, accomplished by using a proper 
clipping function A{-} (e.g., for binary signaling A{-} is the 
sign function). 

For example, linear channels, which appear extensively in 
many applications in communication and data storage systems. 



are characterized by the linear relation 

y = /{x} = Rx + n, 

where n is a if x 1 additive noise vector and R = S^S 
is a positive-definite symmetric matrix, often known as the 
correlation matrix. The N xK matrix S describes the physical 
channel medium while the vector y corresponds to the output 
of a bank of filters matched to the physical channel S. 

Assuming linear channels with AWGN with variance cr^ as 
the ambient noise, the general linear detection rule ([T]l can 
describe known linear detectors. For example [16], [17]: 

• The conventional matched filter (MF) detector is obtained 
by taking A = Ik and b = y. This detector is optimal, 
in the MAP-sense, for the case of zero cross-correlations, 
i.e. , R = Ix, as happens for orthogonal CDMA or when 
there is no ISI effect. 

• The decorrelator (zero forcing equalizer) is achieved by 
substituting A = R and b = y. It is optimal in the 
noiseless case. 

• The linear minimum mean-square error (MMSE) detector 
can also be described by using A = R + (J^Ik- This de- 
tector is known to be optimal when the input distribution 
Px is Gaussian. 

In general, linear detection is suboptimal because of its 
deterministic underlying mechanism {i.e. , solving a given set 
of linear equations), in contrast to other estimation schemes, 
such as MAP or maximum likelihood, that emerge from an 
optimization criterion. 

In [19], linear detection, in its general form ([U, was 
implemented using an efficient message-passing algorithm. 
The linear detection problem was shifted from an algebraic 
to a probabilistic domain. Instead of solving a deterministic 
vector-matrix linear equation, an inference problem is solved 
in a graphical model describing a certain Gaussian distribution 
function. Given the overall channel matrix R and the obser- 
vation vector y, one knows how to write explicitly p(x) and 
the corresponding graph Q with edge potentials ('compatibility 
functions') '4>ij and self-potentials ('evidence') These graph 
potentials are determined according to the following pairwise 
factorization of the Gaussian distribution p(x) 

K 

resulting in ipij{xi,Xj) = exp{—XiRijXj) and 
(l)i{xi) = exp (^biXi — Riixf /2). The set of edges {i,j} 
corresponds to the set of all non-zero entries of A 
for which i > j. Hence, we would like to calculate 
the marginal densities, which must also be Gaussian, 
p{x^) - Af{p, = {'R-^y}^, P^^ = {R~H")' where 
and Pi are the marginal mean and inverse variance (a.k.a. 
precision), respectively. It is shown that the inferred mean ji 
is identical to the desired solution x* ~ R^^y. Table I lists 
the GaBP algorithm update rules. 



TABLE I 

Computing A^^b via GaBP. Online matlab implementation is provided in [25]. 
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v^ompuLe r^ij — /ijj diiu /ijj — 
Set Pki = and fiki = 0, V/c 7^ i. 


2. 


Iterate 


Propagate P^i and fi^i, ^ i such that Aki 7^ 0. 

Compute P,\j = Pu + J2k&m\] and Pr^]{P^^^l^^ + EfcGN(»)\j -Pfe»Mfc«)- 
Compute -AijP^^^Aji and = -P^^^Aijfj.i\j. 


3. 


Check 


If Pij and /ijj did not converge, return to #2. Else, continue to #4. 


4. 


Infer 


= Pa + X]fceN(i) -^fci ' ^ -^i {Piil^ii + SfeGN(i) Pkil^-ki)- 


5. 


Decide 


X., = A{/ii} 



III. Distributed Iterative Computation of the 
MMSE Detector 

In this section, we efficiently extend the applicability of the 
proposed GaBP-based solver for systems with symmetric ma- 
trices [19] to systems with any square (i.e., also nonsymmet- 
ric) or rectangular matrix. We first construct a new symmetric 
data matrix R based on an arbitrary (non-rectangular) matrix 

S e R*-'^" 

^-(s ^ eM('=+")x (2) 

Additionally, we define a new vector of variables x = 
{x^, z^}^ e R(fe+")xi, where x e M'^''^ is the (to be shown) 
solution vector and z e M"^^ is an auxiliary hidden vector, 
and a new observation vector y = {0-^,y-^}-^ G R('^+")'^i. 

Now, we would like to show that solving the symmetric 
linear system Rx = y and taking the first k entries of the 
corresponding solution vector x is equivalent to solving the 
original (not necessarily symmetric) system Rx — y. Note 
that in the new construction the matrix R is sparse again, and 
has only 2nk off-diagonal nonzero elements. When running 
the GaBP algorithm we have only 2nk messages, instead of 
in the previous construction. 

Writing explicitly the symmetric linear system's equations, 
we get 

x + S^z = 0, Sx-^'z = y. 

Thus, 

x = M/"iS^(y-Sx), 

and extracting x we have 

X = (S'^S + *)-iS^y. 

Note, that when the noise level is zero, 5* = Omxm, we get 
the Moore-Penrose pseudoinverse solution 

x-(S^S)-iS^y = Sty. 

IV. Relation to factor graph 

In this section we give an alternate proof of the correctness 
of our construction. Given the inverse covariance matrix R 
defined in (|2]i, and the shift vector x we can derive the 
matching self and edge potentials 



(j)i{xi) = exp(-l/2a;ii?^jXi - Xiyi) 
which is a factorization of the Gaussian system distribution 

i<k i>k ij 

prior on x 

= '[[cxp{-^xf)Y[exp{-^'ii,x'^~x,y,)Y[exp{-x^ S'„ Xj) 

i<k i>k i.j 

Next, we show the relation of our construction to a factor 
graph. We will use a factor graph with k nodes to the left (the 
bits transmitted) and n nodes to the right (the signal received), 
shown in Fig. 1. Using the definition x = {x^, z^}^ € 
]|j(fe+n)xi jjjg vector X represents the k input bits and the 
vector z represents the signal received. Now we can write the 
system probability as: 

p(x) cx / 7V(x; 0, /)A/'(z; Sx, *)dx 

J X 

It is known that the marginal distribution over z is: 

= AA(z;0,S^S + *) 

This distribution is Gaussian, with the following parameters: 

i;(z|x) = (S'^S + *)-iS'^y 

Cov{z\±) ^ (S^S + *)-i 

It is interesting to note that a similar construction was used 
by Frey [26] in his seminal 1999 work when discussing the 
factor analysis learning problem. While his work is beyond 
the scope of this paper, it can be shown that his algorithm can 
be modelled using the GaBP algorithm. 

V. New convergence results 

One of the benefits of using our new construction is that 
we propose a new mechanism to provide future convergence 
results. In the Appendix we prove that Montanari's algorithm 
is an instance of our algorithm, thus our convergence results 
apply to Montanari's algorithm as well. 

We know that if the matrix R is strictly diagonally domi- 
nant, then GaBP converges and the marginal means converge 



CDMA channel linear transformation 




Fig. 1. Factor graph describing the linear channel 



to the true means [27, Claim 4]. Noting that the matrix R is 
symmetric, we can determine the applicabihty of this condition 
by examining its columns. Referring to (4) we see that in 
the first k columns, we have the k CDMA sequences. We 
assume random-spreading binary CDMA sequences which are 
normalized to one. In other words, the absolute sum of each 
column is ^/n. In that case, the matrix R is not diagonally 
dominant (DD). We can add a regularization term of ^/n + e 
to force the matrix R to be DD, but we pay in changing the 
problem. In the next n columns of the matrix R, we have the 
diagonal covariance matrix ^I^ with different noise levels per bit 
in the main diagonal, and zero elsewhere. The absolute sum of 
each column of S is k/^/n, thus when the noise level of each 
bit satisfies '^i > k/y/n, we have a convergence guarantee. 
Note, that the convergence condition is a sufficient condition. 
Based on Montanari's work, we also know that in the large 
system limit, the algorithm converges for binary signaling, 
even in the absence of noise. 

An area of future work is to utilize this observation to 
identify CDMA schemes with matrices S that when fitted into 
the matrix R are either DD, or comply to the spectral radius 
convergence condition of [28]. 

VI. Conclusion 

We presented a novel distributed algorithm for comput- 
ing the MMSE detector for the CDMA multiuser detection 
problem. Our work utilizes the Gaussian Belief Propagation 
algorithm while improving two existing constructions [19], 
[24] in this field. Although we described our algorithm in 
the context of multiuser detection, it has wider applicability. 
For example, it provides an efficient iterative method for 
computing the Moore-Penrose pseudoinverse, and it can also 
be applied to the factor analysis learning problem [26]. 

APPENDIX: MONTANARI'S ALGORITHM IS AN INSTANCE 
OF OUR ALGORITHM 

In this section we show that Montanari's algorithm is an 
instance of our algorithm. Our algorithm is more general. 
First, we allow different noise level for each received bit. 



unlike his work which uses a single fixed noise for the whole 
system. In practice, the bits are transmitted using different 
frequencies, thus suffering from different noise levels. Second, 
the update rules in his paper are fitted only to the randoml- 
spreading CDMA codes, where the matrix A contains only 
values which are drawn uniformly from {—1,1}. Assuming 
binary signalling, he conjectures convergence to the large 
system limit. Our new convergence proof holds for any CDMA 
matrices provided that the absolute sum of the chip sequences 
is one, under weaker conditions on the noise level. Third, 
we propose in [19] an efficient broadcast version for saving 
messages in a broadcast supporting network. 

The probability distribution of the factor graph used by 
Montanari is: 

1^1 ^1 
'^^^v'^ = -:;Wk IT exp(-Trcr^t^^ + jya^^a) n exp(--a;- ) 

• J]^ cxp(----^s„Waa;i)dw 

'i,a 

Extracting the self and edge potentials from the above 
probability distribution: 

^/>,,(x,) = exp(-ix,2) cx AA(x; 0, 1) 

1pia{Xi,UJa) = exp( y=SaiUJayii) OC 7V(x; -y=Sat,0) 

V jV V -'V 

For convenience. Table II provides a translation between the 
notations used in this paper (taken from [19]) and that used 
by Montanari et al. in [24]: 



TABLE II 
Summary of notations 
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At + l) 


precision msg from left to right 


a — '1 


precision msg from right to left 




H — *a 
! a — »-i 


mean msg from left to right 




mean msg from right to left 




Vi 


prior mean of left node 







prior mean of right node 


P 


1 


prior precision of left node 




<72 


prior precision of right node 




Gi 
Li 


posterior mean of node 


Pi 


Li 


posterior precision of node 




•/N 


covariance 






covariance 
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Now we derive Montanari's update rules. We start with the 
precision message from left to right: 

P., 




N \W 

'jSib 1 -jSib 



(0 



N 



By looking at Table 1, it is easy to verify that this precision 
update rule is equivalent to that in #2 of Table I. 

Using the same logic we get the precision message from 
right to left: 




The mean message from left to right is given by 



1 ^ Sg 



N 



-A P-y' 




= -Sb^^a 



The same calculation is done for the mean from right to left: 

-(t+l) Lv ^fca (t) 

A; 

k — m 

Finally, the left nodes calculated the precision and mean by 



g: 



(*+i) 



1 



Sib .(t) 



x,w T. - r--! 



- 1 + TV ''T(*r ' ~ ' i ■ 

The key difference between the two constructions is that 
Montanari uses a directed factor graph while we use an undi- 
rected graphical model. As a consequence, our construction 
provides additional convergence results and simpler update 
rules. 
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